install.packages("rddtools")
install.packages("ggpubr")
install.packages("tidyr")
remove.packages("rlang")
install.packages("rlang")
remove.packages("vctrs")
install.packages("vctrs")
remove.packages("cli")
install.packages("cli")
install.packages("rdd")
install.packages("truncnorm")
library("truncnorm")
library("rdd")
library("rddtools")
library("ggplot2")
library("tidyr")
library("rlang")
library("ggpubr")
library("jtools")
library(stargazer)


rm(list=setdiff(ls(), c("Base", "Base_RDD")))


Base_RDD <- Base %>%
  dplyr::filter(Country=="Brazil",
                ANO == 2010,
                RL != 0) %>%
  group_by(Setor) %>%
  mutate(logRL = log(RL),
         MediaRL_set = mean(logRL),
         DifRL = MediaRL_set - logRL,
         ETR = ETRCOLETADO) %>%
  ungroup()


reg1 <-lm(ETR ~ DifRL, data = Base_RDD)
summary(reg1)

#quanto maior o aumento de receita, maior a alíquota efetiva 

#estimar valores da variável de interesse e testar linearidade

predito <- predict(reg1)
plot(Base_RDD$DifRL, predito, xlab = "Diferencial de Receita", ylab = "Valores previstos de ETR") + abline (v=0)


# A relação entre o diferencial de RL e o ETR não é mais contínua.
# Criou-se uma descontinuidade 


#Gerando o gráfico para visualizar o comportamento da variável y (ETR) em relação à running variable (DifRL)
Graf_ETR_DifRL <- Base_RDD %>% 
  ggplot(aes(x = DifRL, y = ETR)) + 
  geom_point() +
  geom_vline(xintercept = 0.00000, color = "red", size = 1, linetype = "dashed") + 
  annotate("text", x = 0.26415 , y = -0.01819 , label = "Diferença Mínima",
           size=4) +
  labs(y = "ETR",
       x = "DifRL")

print(Graf_ETR_DifRL)

#Gráfico de dispersão
Graf_mesma <- Base_RDD %>%
  dplyr::select(DifRL, ETR) %>%
  mutate(limite = as.factor(ifelse(DifRL >= 0.00000, 1, 0))) %>%
  ggplot(aes(x = DifRL, y = ETR)) +
  geom_point(aes(color = limite)) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_brewer(palette = "Accent") +
  guides(color = FALSE) +
  geom_vline(xintercept = 0.00000, color = "red",
             size = 1, linetype = "dashed") +
  labs(y = "ETR",
       x = "DifRL")

plot(Graf_mesma)

#Identificando a amplitude da amostra
summary(Base_RDD$DifRL)
summary(Base_RDD$ETR)

#Estimando o efeito usando a função lm
Base_RDD$limite <- ifelse(Base_RDD$DifRL >= 0.00000, 1, 0)
Mod_RDD_lm <- lm(ETR ~ I(DifRL - 0.00000) + limite, data = Base_RDD)

summary(Mod_RDD_lm) 


# Em resumo, com base nos resultados apresentados, não há evidências 
# de uma diferença entre o grupo tratamento e o grupo controle

#Normalização e estimação

Base_RDD$DifRL_left <- ifelse(Base_RDD$DifRL <=0, Base_RDD$DifRL-0, 0)
Base_RDD$DifRL_right <- ifelse(Base_RDD$DifRL >0, Base_RDD$DifRL-0, 0)
Base_RDD$limite <- ifelse(Base_RDD$DifRL >= 0.00000, 1, 0)
reg2 <- lm(ETR ~ DifRL_left + DifRL_right + limite, Base_RDD)
summary(reg2)
predito2 <- predict(reg2)
plot(Base_RDD$DifRL, predito2, xlab = "Diferencial de Receita", ylab = "Valores previstos de ETR") + abline (v=0)


#Estimando usando o pacote RDD Tools
Mod_RDD_Tools <- rdd_data(y = Base_RDD$ETR, 
         x = Base_RDD$DifRL, 
         cutpoint = 0.00000)

Mod_RDD_Tools2 <- rdd_reg_lm(Mod_RDD_Tools, slope = "same")

summary(Mod_RDD_Tools2)

# slope: A opção "same" indica  a inclinação (slope) 
# da reta de regressão com a mesma inclinação antes e depois do ponto de corte.


#Estimando o efeito com diferentes inclinações usando a função lm
Mod_RDD_Dif_lm <- Base_RDD %>% 
  mutate(limite = ifelse(DifRL >= 0.00000, 1, 0)) 
Mod_RDD_Dif_lm <- lm(ETR ~ limite + I(DifRL - 0.00000)+ limite:I(DifRL - 0.00000), data = Base_RDD)

summary(Mod_RDD_Dif_lm)

#Estimando usando o pacote RDD Tools
Mod_RDD_Dif_Tools <- rdd_data(y = Base_RDD$ETR, 
         x = Base_RDD$DifRL, 
         cutpoint = 0.00000)

Mod_RDD_Dif_Tools2 <- rdd_reg_lm(Mod_RDD_Dif_Tools, slope = "separate")

# slope = "separate": indica que a análise está 
# considerando inclinações separadas em cada lado do ponto de corte. 

summary(Mod_RDD_Dif_Tools2)
print(Mod_RDD_Dif_Tools2)

#Gráfico de dispersão para diferentes inclinações

Graf_Dif <- Base_RDD %>%
  dplyr::select(DifRL, ETR) %>%
  mutate(limite = as.factor(ifelse(DifRL >= 0.00000, 1, 0))) %>%
  ggplot(aes(x = DifRL, y = ETR, color = limite)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_brewer(palette = "Accent") +
  guides(color = FALSE) +
  geom_vline(xintercept = 0.00000, color = "red",
             size = 1, linetype = "dashed") +
  labs(y = "ETR",
       x = "DifRL")

ggarrange(Graf_ETR_DifRL, Graf_mesma, Graf_Dif)

summary(Base_RDD$DifRL, details)

Graf_ETR_DifRL

#Modelo de sensibilidade

Mod_Sensibilidade <- Base_RDD %>%
  dplyr::filter(DifRL >= -1 & DifRL <= 1) %>%
  mutate(limite = ifelse(DifRL >= 0.00000, 1, 0))

Mod_Sensibilidade <- lm(ETR ~ limite + I(DifRL - 0.00000) + limite:I(DifRL - 0.00000),data = Base_RDD)

summary(Mod_Sensibilidade)


#Gráfico de dispersão da análise de sensibilidade
Graf_Sensibilidade <- Base_RDD %>%
  dplyr::filter(DifRL >= -1 & DifRL <= 1) %>%
  dplyr::select(DifRL, ETR) %>%
  mutate(limite = as.factor(ifelse(DifRL >= 0.00000, 1, 0))) %>%
  ggplot(aes(x = DifRL, y = ETR, color = limite)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_brewer(palette = "Accent") +
  guides(color = FALSE) +
  geom_vline(xintercept = 0.00000, color = "red",
             size = 1, linetype = "dashed") +
  labs(y = "ETR",
       x = "DifRL")


#Alterando a forma funcional
Mod_Sens_Quad <- Base_RDD %>%
  mutate(limite = ifelse(DifRL >= 0.00000, 1, 0))

  
Mod_Sens_Quad <- lm(ETR ~ limite + I(DifRL - 0.00000)  + I((DifRL - 0.00000)^2) + limite:I(DifRL - 0.00000) +
       limite:I((DifRL - 0.00000)^2), data = Base_RDD)

summary(Mod_Sens_Quad)

#Gráfico de dispersão
Graf_Quad <- Base_RDD %>%
  dplyr::select(DifRL, ETR) %>%
  mutate(limite = as.factor(ifelse(DifRL >= 0.00000, 1, 0))) %>%
  ggplot(aes(x = DifRL, y = ETR, color = limite)) +
  geom_point() +
  geom_smooth(method = "lm",
              formula = y ~ x + I(x ^ 2),
              se = FALSE) +
  scale_color_brewer(palette = "Accent") +
  guides(color = FALSE) +
  geom_vline(xintercept = 0.00000, color = "red",
             size = 1, linetype = "dashed") +
  labs(y = "ETR",
       x = "DifRL")


export_summs(Mod_RDD_lm, Mod_RDD_Dif_lm, Mod_Sensibilidade,Mod_Sens_Quad)
ggarrange(Graf_ETR_DifRL, Graf_mesma, Graf_Dif, Graf_Sensibilidade, Graf_Quad)


#### FUZZY RDD

df <-data.frame(A = rtruncnorm(10000, a = 0, mean = 1000, sd = 500))

df <- df %>%
  mutate(prob_tratamento = ifelse(A<500,0.8,0),
         tratamento = rbinom(10000, 1, prob = prob_tratamento),
         y = A + 150*tratamento + rnorm(10000, 100, 50),
         tratamento = as.factor(tratamento))

ggplot(df, aes(x= A, y= y)) + 
  geom_point(aes(color = tratamento))

fuzzy_rdd = rdd_data(y=y,
                     x=A,
                     z = df$prob_tratamento,
                     cutpoint = 500, data = df) %>%
  rdd_reg_lm(slope = "separate")
summary(fuzzy_rdd)
    
